Resource selection function



All data aggregated but two different fundamental assumption sets depending on how constraints associated with unused points. Unused can be fully randomly sampled, i.e. could in theory overlap with potential animal relocation at a different point in space or time. Alternatively, unused points sampled from environment can be constrainted to ensure that they never overlap with another relocation instance.


Models

#data with replace means can overlap (same assumption as work in Carrizo)
#data <- read_csv("data/data_with_replace.csv")
data <- read_csv("data/tidy_data.csv") 
data <- data %>%
  select(rep, status, year, ground, mesohabitat, Aspect, Slope, Elev, KDF, Shrub_Cov, NDVI, Solar, aspect.cat)
data$year <-  as.character(data$year)
data
## # A tibble: 6,708 x 13
##      rep status year  ground mesohabitat Aspect Slope  Elev     KDF
##    <dbl>  <dbl> <chr> <chr>  <chr>        <dbl> <dbl> <dbl>   <dbl>
##  1     2      0 2016  above  open         217.  3.58    708 3.24e-3
##  2     3      0 2018  below  open         218.  4.08    709 3.77e-3
##  3     7      0 2018  above  shrub        247.  2.73    697 7.69e-5
##  4     8      0 2018  above  shrub        202.  3.85    710 3.72e-3
##  5     9      0 2018  above  shrub        225   0.506   698 3.26e-5
##  6    11      0 2016  below  open         225   3.54    709 3.62e-3
##  7    16      0 2018  above  shrub        288.  1.13    692 6.90e-6
##  8    21      0 2016  above  open         217.  3.58    708 3.24e-3
##  9    24      0 2017  above  shrub         63.4 1.60    698 4.10e-6
## 10    27      0 2018  above  shrub        297.  3.20    698 0.     
## # ... with 6,698 more rows, and 4 more variables: Shrub_Cov <dbl>,
## #   NDVI <dbl>, Solar <dbl>, aspect.cat <chr>
#function
#global availability (m=0) and bootstrap (B=99)
# m = 0 indicates sample from anywhere in dataframe versus structure or blocks
model.rspf <- function(x) {rspf(status ~ Shrub_Cov + NDVI + Slope + Elev + Aspect + Solar, x, m=0, B = 99, model = TRUE)
}

m7 <- model.rspf(data)
summary(m7)
## 
## Call:
## rspf(formula = status ~ Shrub_Cov + NDVI + Slope + Elev + Aspect + 
##     Solar, data = x, m = 0, B = 99, model = TRUE)
## 
## Resource Selection Probability Function (Logistic RSPF) model
## Non-matched Used-Available design
## Maximum Likelihood estimates
## with Nonparametric Bootstrap standard errors (B = 99)
## 
## Fitted probabilities:
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.8263  0.9963  0.8236  0.9997  1.0000 
## 
## Coefficients (logit link):
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  2.228e+01  5.510e-01  40.440  < 2e-16 ***
## Shrub_Cov    9.580e+00  3.521e-01  27.205  < 2e-16 ***
## NDVI         9.132e+00  7.249e-01  12.597  < 2e-16 ***
## Slope        2.428e-01  3.575e-02   6.792  1.1e-11 ***
## Elev        -2.663e-01  1.322e-02 -20.138  < 2e-16 ***
## Aspect      -3.043e-04  1.666e-03  -0.183    0.855    
## Solar        4.469e-04  2.564e-05  17.429  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Log-likelihood: -2.647e+04 
## BIC = 5.3e+04 
## 
## Hosmer and Lemeshow goodness of fit (GOF) test:
## X-squared = 292.4, df = 8, p-value < 2.2e-16
plot(m7)

mep(m7) # marginal effects similar to plot but with CIs

kdepairs(m7) # 2D kernel density estimates

#split-map to explore key levels####
#year####
m.year <- data %>% 
  split(.$year) %>% 
  map(~ model.rspf(.))
m.year
## $`2016`
## 
## Call:
## rspf(formula = status ~ Shrub_Cov + NDVI + Slope + Elev + Aspect + 
##     Solar, data = x, m = 0, B = 99, model = TRUE)
## 
## Resource Selection Probability Function (Logistic RSPF) model
## Non-matched Used-Available design
## Maximum Likelihood estimates
## 
## Coefficients (logit link):
## (Intercept)    Shrub_Cov         NDVI        Slope         Elev  
##  -2.059e+00    1.653e+01    4.346e+00   -6.792e-02    1.188e-02  
##      Aspect        Solar  
##  -1.447e-03   -2.424e-05  
## 
## 
## $`2017`
## 
## Call:
## rspf(formula = status ~ Shrub_Cov + NDVI + Slope + Elev + Aspect + 
##     Solar, data = x, m = 0, B = 99, model = TRUE)
## 
## Resource Selection Probability Function (Logistic RSPF) model
## Non-matched Used-Available design
## Maximum Likelihood estimates
## 
## Coefficients (logit link):
## (Intercept)    Shrub_Cov         NDVI        Slope         Elev  
##  -5.921e+01    1.347e+00    2.099e+01   -8.013e-02   -6.508e-03  
##      Aspect        Solar  
##   1.274e-03    1.497e-04  
## 
## 
## $`2018`
## 
## Call:
## rspf(formula = status ~ Shrub_Cov + NDVI + Slope + Elev + Aspect + 
##     Solar, data = x, m = 0, B = 99, model = TRUE)
## 
## Resource Selection Probability Function (Logistic RSPF) model
## Non-matched Used-Available design
## Maximum Likelihood estimates
## 
## Coefficients (logit link):
## (Intercept)    Shrub_Cov         NDVI        Slope         Elev  
##   6.967e+01    2.933e+00    1.399e+01    7.802e-02   -2.341e-01  
##      Aspect        Solar  
##   3.633e-05    2.521e-04
#ground####
m.ground <- data %>% 
  split(.$ground) %>% 
  map(~ model.rspf(.))
m.ground
## $above
## 
## Call:
## rspf(formula = status ~ Shrub_Cov + NDVI + Slope + Elev + Aspect + 
##     Solar, data = x, m = 0, B = 99, model = TRUE)
## 
## Resource Selection Probability Function (Logistic RSPF) model
## Non-matched Used-Available design
## Maximum Likelihood estimates
## 
## Coefficients (logit link):
## (Intercept)    Shrub_Cov         NDVI        Slope         Elev  
##  10.2942123   10.2087814    5.4884795    0.3357263   -0.3496654  
##      Aspect        Solar  
##   0.0005631    0.0006408  
## 
## 
## $below
## 
## Call:
## rspf(formula = status ~ Shrub_Cov + NDVI + Slope + Elev + Aspect + 
##     Solar, data = x, m = 0, B = 99, model = TRUE)
## 
## Resource Selection Probability Function (Logistic RSPF) model
## Non-matched Used-Available design
## Maximum Likelihood estimates
## 
## Coefficients (logit link):
## (Intercept)    Shrub_Cov         NDVI        Slope         Elev  
##  24.6215758    8.3778257   10.8635911    0.1259252   -0.2341786  
##      Aspect        Solar  
##  -0.0003167    0.0003786
#mesohabitat
m.meso <- data %>% 
  split(.$mesohabitat) %>% 
  map(~ model.rspf(.))
m.meso
## $open
## 
## Call:
## rspf(formula = status ~ Shrub_Cov + NDVI + Slope + Elev + Aspect + 
##     Solar, data = x, m = 0, B = 99, model = TRUE)
## 
## Resource Selection Probability Function (Logistic RSPF) model
## Non-matched Used-Available design
## Maximum Likelihood estimates
## 
## Coefficients (logit link):
## (Intercept)    Shrub_Cov         NDVI        Slope         Elev  
##  36.6829900   -0.4861748    3.6675175    0.2149919   -0.2715408  
##      Aspect        Solar  
##   0.0018867    0.0004218  
## 
## 
## $shrub
## 
## Call:
## rspf(formula = status ~ Shrub_Cov + NDVI + Slope + Elev + Aspect + 
##     Solar, data = x, m = 0, B = 99, model = TRUE)
## 
## Resource Selection Probability Function (Logistic RSPF) model
## Non-matched Used-Available design
## Maximum Likelihood estimates
## 
## Coefficients (logit link):
## (Intercept)    Shrub_Cov         NDVI        Slope         Elev  
##  -2.843e+00    2.377e+01    8.713e+00   -8.904e-02   -1.748e-02  
##      Aspect        Solar  
##   2.074e-03    2.966e-05

Effect size estimate

A few instances in literature calculate log odds ratio to estimate effect sizes. Contrast of used by unused probabilities.

#build dataframe
#model <- as_tibble(m7$model)
#fitted_values <- as_tibble(m7$fitted.values)
#odds <- bind_cols(model, fitted_values)
#odds
#write_csv(odds, "data/odds.csv")
#library(esc)
#effect_sizes(odds, t = value, fun = "esc_bin_prop", es.type = "or")

#just to run this section only
library(tidyverse)
#library(scales)
#odds <- read_csv("data/odds.csv")

#log odds ratios####
#log(p(A),p(not-A)) typically for A as event or cell frequency
#log((p1/(1-p1))/(p2/(1-p2))) for proportion successful in each group

odds_df <- read_csv("data/log_odds.csv") %>%
  select(-aspect)

odds_df_long <- odds_df %>%
  gather(factor, measure, 1:5)

ggplot(odds_df_long, aes(measure, log_odds)) +
  geom_smooth(method = lm, color = "darkgreen") +
  facet_wrap(~factor, scales = "free") +
  labs(x = "", y = "log odds ratio")

library(broom)
m <- odds_df_long %>%
  group_by(factor) %>%
  do(fit = lm(log_odds~measure, data =.)) 
tidy(m, fit)
## # A tibble: 10 x 6
## # Groups:   factor [5]
##    factor      term           estimate  std.error statistic  p.value
##    <chr>       <chr>             <dbl>      <dbl>     <dbl>    <dbl>
##  1 elevation   (Intercept)   78.9       1.43          55.4  0.      
##  2 elevation   measure       -0.104     0.00200      -51.8  0.      
##  3 NDVI        (Intercept)    4.64      0.388         12.0  2.82e-32
##  4 NDVI        measure        1.87      1.70           1.10 2.72e- 1
##  5 shrub_cover (Intercept)    4.94      0.0517        95.5  0.      
##  6 shrub_cover measure        5.38      1.54           3.50 4.68e- 4
##  7 slope       (Intercept)    5.27      0.0760        69.3  0.      
##  8 slope       measure       -0.0672    0.0223        -3.02 2.54e- 3
##  9 solar       (Intercept)  165.       10.9           15.2  1.84e-50
## 10 solar       measure       -0.000422  0.0000286    -14.7  1.29e-47

Interpretations

  1. No difference between overlapping and non-overlapping unused sample points (EDA).
  2. Best models include all key environnmental drivers.
  3. Shrub cover, NDVI, and solar all have positive effects on probability of sampling animals.
  4. Aspect is not statistically significant (and is redundant with solar exposure effects).
  5. Increasing elevation decreasses likelihood of sampling animals.
  6. The relative importance of different drivers can vary annually but were consistent in general effects.
  7. Shrub cover predicted probabilities above and below ground.
  8. Shrub cover (downscaled data) as a predictor was more important for shrub relocations confirming data derivation process.